#include <iostream>
#include <vector>
#include <ctime>
#include <numeric>

#include "tensor.h"
#include "covariance.h"
#include "mtrand.h"
#include "dataset.h"
#include "pdf.h"
#include <vector>
#include <sstream>

#include "sdataset.h"
#include "spdf.h"

#include "TFile.h"
#include "TChain.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TProfile2D.h"
#include "TF1.h"
#include "TF2.h"
#include "TF3.h"
#include "TLorentzVector.h"
#include "TGraphErrors.h"
#include "TGraph2DErrors.h"
#include "make_list.h"

using namespace std;

//initialize random number generator
MTRand_int32 irand(time(0));

const double PI=3.14159265359;
//const double exp10=-2.30258509299;

//function to get the correlation from MC
double MC_correlation(const vector<PDF::MC_Event>& input)
{
  assert(input.size()>0);

  //output should be two variables
  assert(input[0].value.size() == 2);
  
  double mean_xy=0;
  double mean_x=0;
  double mean_y=0;
  double sum_weight=0;

  for(unsigned int j=0; j<input.size(); ++j)
    {
      double temp_x=input[j].value[0];
      double temp_y=input[j].value[1];
      double temp_weight=input[j].weight;
      sum_weight+=temp_weight;
      mean_xy+=temp_x*temp_y*temp_weight;
      mean_x+=temp_x*temp_weight;
      mean_y+=temp_y*temp_weight;
    }

  //sum of weight should be greater than zero
  assert(sum_weight > 0);

  sum_weight=1.0/sum_weight;
  
  return (mean_xy*sum_weight - mean_x*mean_y*sum_weight*sum_weight);
}

double MC_correlation_norm(const vector<PDF::MC_Event>& input)
{
  assert(input.size()>0);

  //output should be two variables
  assert(input[0].value.size() == 2);
  
  double mean_xy=0;
  double mean_xx=0;
  double mean_yy=0;

  double mean_x=0;
  double mean_y=0;
  double sum_weight=0;

  for(unsigned int j=0; j<input.size(); ++j)
    {
      double temp_x=input[j].value[0];
      double temp_y=input[j].value[1];
      double temp_weight=input[j].weight;
      sum_weight+=temp_weight;
      mean_xy+=temp_x*temp_y*temp_weight;
      mean_xx+=temp_x*temp_x*temp_weight;
      mean_yy+=temp_y*temp_y*temp_weight;
      mean_x+=temp_x*temp_weight;
      mean_y+=temp_y*temp_weight;
    }

  //sum of weight should be greater than zero
  assert(sum_weight > 0);
  sum_weight=1.0/sum_weight;

  double std_xx=sqrt(mean_xx*sum_weight - mean_x*sum_weight*mean_x*sum_weight);
  double std_yy=sqrt(mean_yy*sum_weight - mean_y*sum_weight*mean_y*sum_weight);
  
  assert(std_xx>0);
  assert(std_yy>0);
  
  return
    (mean_xy*sum_weight - mean_x*mean_y*sum_weight*sum_weight)/
    (std_xx*std_yy);
}

//compute the reweight
double MC_reweight(PDF::MC_Event& evt1,
		   PDF::MC_Event& evt2, 
		   double correlation, int nMC)
{
  assert(evt1.value.size()>0);
  assert(evt2.value.size()>0);
  assert(evt1.GetExtra().size()>0);
  assert(evt2.GetExtra().size()>0);
  
  return nMC*evt1.weight*nMC*evt2.weight +
    correlation*evt1.GetExtra()[0]*evt2.GetExtra()[0];
}

//compute rescaling factor
double MC_rescale(vector<PDF::MC_Event>& input1, 
		  vector<PDF::MC_Event>& input2, 
		  double correlation, double& rescale_reweight)
{
  assert(input1.size() == input2.size());	 
  double sum_weight=0;
  double sum_reweight=0;

  for(unsigned int j=0; j<input1.size(); ++j)
    {
      sum_weight+=input1[j].weight*input2[j].weight;
      
      sum_reweight+= 
	MC_reweight(input1[j], input2[j],
		    correlation, input1.size());
    }

  assert(sum_weight > 0);
  assert(sum_reweight != 0);
  sum_weight=1.0/sum_weight;
  rescale_reweight=1.0/sum_reweight;
  return sum_weight;
}



float DeltaR(float eta, float phi)
{
  if(fabs(phi)> PI)
    phi= PI-fabs(phi);
  return sqrt(eta*eta + phi*phi);
}

int main(int argc, char **argv)
{

  //flag = 2jet, 4jet, 3jet, 2jet-exclusive, 3jet-exclusive, 4jet-exclusive"
  //ptcut = cuts to define jets
  //default is 2jet

  if(argc < 4+1)
    {
      cout<<"Usage: ./example_sherpa input.root output.root nevt nproc "
	  <<"[flag] [ptcut]"<<endl;
      throw;
    }
  
  //option to change the number of jets
  int njet=2;
  bool b_exclusive=false;
  int ptcut=100;
  if(argc>4+1)
    {
      string snjet(argv[5]);
      if(snjet.find("3jet")!=string::npos)
	njet=3;
      else if(snjet.find("4jet")!=string::npos)
	njet=4;

      cout<<njet<<" jets";
      if(snjet.find("exclusive")!=string::npos)
	{
	  b_exclusive=true;
	  cout<<" exclusive sample";
	}

      if(argc>5+1)
	{
	  ptcut= atoi(argv[6]);
	}
      
      cout<<" with "<<ptcut<<" GeV pt cut"<<endl;
    }
  
  cout<<"Implicit input from PDF_* files"<<endl;
  
  int nevt=atoi(argv[3]);
  int nproc=atoi(argv[4]);

  assert(nevt>0);
  assert(nproc>=0);

  //open rootfile
  TChain chain("Ntuple");
  chain.Add(argv[1]);

  if(chain.GetEntries()<nevt*nproc)
    {
      cout<<"End of file, exit..."<<endl;
      return 1;
    }

  //open output file
  TFile outfile(argv[2], "RECREATE");

  //declare variables
  const int maxsize=10;
  Int_t           size;
  Double_t         w;   //[size]
  Float_t         pt[maxsize];   //[size]
  Float_t         eta[maxsize];
  Float_t         phi[maxsize];
  Float_t         m[maxsize];   //[size]
  Float_t         t1[maxsize];
  Float_t         t2[maxsize];
  Float_t         t3[maxsize];


  //read the input file
  chain.SetBranchAddress("size", &size);
  chain.SetBranchAddress("w", &w);
  chain.SetBranchAddress("pt", pt);
  chain.SetBranchAddress("eta", eta);
  chain.SetBranchAddress("phi", phi);
  chain.SetBranchAddress("m", m);
  chain.SetBranchAddress("tau1", t1);
  chain.SetBranchAddress("tau2", t2);
  chain.SetBranchAddress("tau3", t3);

  //declare some histograms
  map<string, TH1*> hist;
  hist["t21_1"]=new TH1D("t21_1", "Leading Jet tau 21", 20, 0, 1.0);
  hist["t21_2"]=new TH1D("t21_2", "Second Jet tau 21", 20, 0, 1.0);

  hist["t21_1_MC"]=new TH1D("t21_1_MC", "Leading Jet tau 21 (KDE)", 20, 0, 1.0);
  hist["t21_1_MCerr"]=new TH1D("t21_1_MCerr", "Leading Jet tau 21 (KDE)",
			       20, 0, 1.0);

  hist["t21_2_MC"]=new TH1D("t21_2_MC", "Second Jet tau 21 (KDE)", 20, 0, 1.0);
  hist["t21_2_MCerr"]=new TH1D("t21_2_MCerr", "Second Jet tau 21 (KDE)",
			       20, 0, 1.0);

  hist["T21"]=new TH1D("T21", 
		       "T_{21}", 20, 0, 1.0);

  hist["T21_MC"]=new TH1D("T21_MC", 
			  "T_{21} from KDE", 20, 0, 1.0);
  hist["T21_MCerr"]=new TH1D("T21_MC_err", 
			     "T_{21} err from KDE", 20, 0, 1.0);

  hist["T21_MC_corr"]=new TH1D("T21_corr", 
				  "T_{21} from KDE (corrected)", 20, 0, 1.0);
  hist["T21_MCerr_corr"]=new TH1D("T21_corr_err", 
				  "T_{21} err from KDE (corrected)", 20, 0, 1.0);

  hist["m1"]=new TH1D("m1", "mass for leading jet", 40, 50, 500);
  hist["m2"]=new TH1D("m2", "mass for second jet", 40, 50, 500);

  hist["MJ"]=new TH1D("MJ", "M_{J}", 40, 50, 500);
  hist["MJ_4jet"]=new TH1D("MJ_4jet", "4jet M_{J}", 40, 50, 1000);


  hist["MJ_4jet_MC"]=new TH1D("MJ_4jet_MC", "4jet M_{J} (KDE)", 40, 50, 1000);
  hist["MJ_4jet_MCerr"]=new TH1D("MJ_4jet_MCerr", "4jet M_{J} (KDE)", 40, 50, 1000);

  hist["MJ_MC"]=new TH1D("MJ_MC", "M_{J} (KDE)", 40, 50, 500);
  hist["MJ_MCerr"]=new TH1D("MJ_MCerr", "M_{J} (KDE)", 40, 50, 500);

  hist["MJ_2jet_MC"]=new TH1D("MJ_2jet_MC", "M_{J} (KDE)", 40, 50, 500);
  hist["MJ_2jet_MCerr"]=new TH1D("MJ_2jet_MCerr", "M_{J} (KDE)", 40, 50, 500);
  
  hist["MJ_MC_corr"]=new TH1D("MJ_corr", "M_{J} (KDE)", 40, 50, 500);
  hist["MJ_MCerr_corr"]=new TH1D("MJ_corr_err", "M_{J} (KDE)", 40, 50, 500);

  hist["MJ_2jet_MC_corr"]=new TH1D("MJ_2jet_corr", "M_{J} (KDE)", 40, 50, 500);
  hist["MJ_2jet_MCerr_corr"]=new TH1D("MJ_2jet_corr_err", "M_{J} (KDE)", 40, 50, 500);

  hist["m1_MC"]=new TH1D("m1_MC", "mass for leading jet (KDE)", 40, 50, 500);
  hist["m2_MC"]=new TH1D("m2_MC", "mass for second jet (KDE)", 40, 50, 500);
  
  hist["m1_MCerr"]=new TH1D("m1_MCerr", "mass for leading jet (KDE)", 40, 50, 500);
  hist["m2_MCerr"]=new TH1D("m2_MerrC", "mass for second jet (KDE)", 40, 50, 500);

  
  //write all histograms
  map<string, TH1*>::iterator it;
  for(it=hist.begin(); it!=hist.end(); ++it)
    {
      it->second->Sumw2();
    }

  
  //here compute come PDF
  PDF PDF_tau21_2jet_dR;
  if(!PDF_tau21_2jet_dR.Read("PDF_tau21_2jet_dR"))
    {
      cout<<"PDF READING ERROR"<<endl;
      return 0;
    }

  PDF PDF_tau21_2jet_mass;
  if(!PDF_tau21_2jet_mass.Read("PDF_tau21_2jet_mass"))
    {
      cout<<"PDF READING ERROR"<<endl;
      return 0;
    }

  PDF PDF_mpt_2jet_mass;
  if(!PDF_mpt_2jet_mass.Read("PDF_mpt_2jet_mass"))
    {
      cout<<"PDF READING ERROR"<<endl;
      return 0;
    }

  PDF PDF_mpt_2jet_dR;
  if(!PDF_mpt_2jet_dR.Read("PDF_mpt_2jet_dR"))
    {
      cout<<"PDF READING ERROR"<<endl;
      return 0;
    }


  PDF PDF_mpt_2jet_pt1_dR;
  if(!PDF_mpt_2jet_pt1_dR.Read("PDF_mpt_2jet_pt1_dR"))
    {
      cout<<"PDF READING ERROR"<<endl;
      return 0;
    }

  PDF PDF_mpt12_pt12;
  if(!PDF_mpt12_pt12.Read("PDF_mpt12_pt12"))
    {
      cout<<"PDF READING ERROR"<<endl;
      return 0;
    }
  
  PDF PDF_tau21_pt;
  if(!PDF_tau21_pt.Read("PDF_tau21_pt"))
    {
      cout<<"PDF READING ERROR"<<endl;
      return 0;
    }

  PDF PDF_mpt_pt;
  if(!PDF_mpt_pt.Read("PDF_mpt_pt"))
    {
      cout<<"PDF READING ERROR"<<endl;
      return 0;
    }


  //get MC_Generator for computation
  PDF::MC_Generator MC_tau21_2jet_dR= PDF_tau21_2jet_dR.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT, PDF::INPUT));

  PDF::MC_Generator MC_mpt_2jet_pt1_dR= PDF_mpt_2jet_pt1_dR.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT, PDF::INPUT, PDF::INPUT));

  PDF::MC_Generator MC_tau21_2jet_mass= PDF_tau21_2jet_mass.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT, PDF::INPUT));

  PDF::MC_Generator MC_mpt_2jet_mass= PDF_mpt_2jet_mass.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT, PDF::INPUT));

  PDF::MC_Generator MC_mpt_2jet_dR= PDF_mpt_2jet_dR.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT, PDF::INPUT));

  PDF::MC_Generator MC_tau21_pt= PDF_tau21_pt.Generator
    (make_list<PDF::MC_Flag>(PDF::OUTPUT, PDF::INPUT));

  PDF::MC_Generator MC_mpt_pt= PDF_mpt_pt.Generator
    (make_list<PDF::MC_Flag>(PDF::OUTPUT, PDF::INPUT));

  PDF::MC_Generator MC_mpt12_pt12= PDF_mpt12_pt12.Generator
    (make_list<PDF::MC_Flag>(PDF::OUTPUT, PDF::OUTPUT,
			     PDF::INPUT, PDF::INPUT));
  
  //double weight_cut=total_weight*0.05/chain.GetEntries();
  //double weight_cut_high=total_weight*0.05;
  //loop over events a second time 

  double low_weight=0;
  double total_weight=0;
  int start=nevt*nproc;
  int end=(nproc+1)*nevt;

  

  DataSet DataSet_mpt_2jet(make_list<double>(0, 0),
			   make_list<double>(2.0,2.0),
			   make_list<double>(0.02,0.02));

  DataSet DataSet_mpt_2jet_KDE(make_list<double>(0, 0),
			       make_list<double>(2.0,2.0),
			       make_list<double>(0.02,0.02));

  DataSet DataSet_mpt_2jet_KDE_corr(make_list<double>(0, 0),
				    make_list<double>(2.0,2.0),
				    make_list<double>(0.02,0.02));

  DataSet DataSet_tau21_2jet(make_list<double>(0, 0),
			     make_list<double>(1.0,1.0),
			     make_list<double>(0.01,0.01));
  
  DataSet DataSet_tau21_2jet_KDE(make_list<double>(0, 0),
				 make_list<double>(1.0,1.0),
				 make_list<double>(0.01,0.01));

  DataSet DataSet_tau21_2jet_KDE_corr(make_list<double>(0, 0),
				      make_list<double>(1.0,1.0),
				      make_list<double>(0.01,0.01));

  bool b_1weight=false;

  for(Long64_t i=start; i<chain.GetEntries() && i<end; i++)
    {
      if((i-start)%10000==0)
	cout<<"\r\033[K"<<((i-start)*10000/nevt)*0.01<<"% computed..."
	    <<flush;

      chain.GetEntry(i);

      /*
	if(w<weight_cut)
	{
	  low_weight+=w;
	  continue;
	}
	if(w>weight_cut_high) continue;
      */

      if(b_1weight)
	w=1.0;

      //njet cut
      if(size<njet) continue;
      
      bool b_cut=true;

      for(int j=0; j<size && j<njet; ++j)
	{
	  //some sanity cuts
	  if(t1[j] < 1e-5 || t2[j] < 1e-5 || pt[j] < ptcut || m[j]<20) 
	    {
	      b_cut=false;
	      break;
	    }
	}

      if(b_exclusive)
	for(int j=njet-1; j<size; ++j)
	    {
	      if( !(t1[j] < 1e-5 || t2[j] < 1e-5 || m[0]<20) &&
		  pt[j] > ptcut)
		b_cut=false;
	    }

      //event does not pass njet criterion
      if(!b_cut) continue;

      if(njet>=4)
	hist["MJ_4jet"]->Fill(m[0]+m[1]+m[2]+m[3], w);
      
      double dR=DeltaR(eta[0]-eta[1], phi[0]-phi[1]);
      TLorentzVector j1,j2;
      j1.SetPtEtaPhiM(pt[0],eta[0],phi[0],m[0]);
      j2.SetPtEtaPhiM(pt[1],eta[1],phi[1],m[1]);
      double mjj=(j1+j2).M();

      double t21_1=t2[0]/t1[0];
      double t21_2=t2[1]/t1[1];
      double mpt_1=-log10(m[0]/pt[0]);
      double mpt_2=-log10(m[1]/pt[1]);

      DataSet_mpt_2jet.Fill(make_list<double>(mpt_1, mpt_2), w);
      DataSet_tau21_2jet.Fill(make_list<double>(t21_1, t21_2), w);


      hist["t21_1"]->Fill(t21_1, w);
      hist["t21_2"]->Fill(t21_2, w);

      hist["T21"]->Fill(sqrt(t21_1*t21_2), w);


      hist["MJ"]->Fill(m[0]+m[1], w);
      hist["m1"]->Fill(m[0], w);
      hist["m2"]->Fill(m[1], w);

      //generate tau21 correlation
      vector<PDF::MC_Event> MC_corr=
	MC_tau21_2jet_mass.GenerateFull(make_list<double>(log(mjj/50.)));
      
      double correlation=MC_correlation(MC_corr);
      
      

      int nGenerate=500;

      //generate some toy MC
      vector<PDF::MC_Event> MC_jet1=
	MC_tau21_pt.GenerateGrad(make_list<double>(log(pt[0]/50.)), nGenerate);

      vector<PDF::MC_Event> MC_jet2=
	MC_tau21_pt.GenerateGrad(make_list<double>(log(pt[1]/50.)), nGenerate);
      
      //since we are doing something complicated
      //we need to rescale all weight to 1
      double rescale=0;
      double rescale_reweight=0;

      //compute rescaling
      rescale=MC_rescale(MC_jet1, MC_jet2, correlation, rescale_reweight);

      for(unsigned int j=0; j<MC_jet1.size(); ++j)
	{
	  hist["t21_1_MC"]->Fill(MC_jet1[j].value[0], w*MC_jet1[j].weight);
	  hist["t21_1_MCerr"]->Fill(MC_jet1[j].value[0], w*MC_jet1[j].err);
	  hist["t21_2_MC"]->Fill(MC_jet2[j].value[0], w*MC_jet2[j].weight);
	  hist["t21_2_MCerr"]->Fill(MC_jet2[j].value[0], w*MC_jet2[j].err);
	}

      //now loop over the MC_event
      for(unsigned int j=0; j<MC_jet1.size(); ++j)
	{
	  double temp_t21_1=MC_jet1[j].value[0];
	  double temp_t21_2=MC_jet2[j].value[0];
	  double temp_T21=sqrt(temp_t21_1*temp_t21_2);

	  double temp_weight=MC_jet1[j].weight*MC_jet2[j].weight*rescale;
	  if(temp_weight==0) continue;
	  double temp_err1=MC_jet1[j].err/MC_jet1[j].weight;
	  double temp_err2=MC_jet2[j].err/MC_jet2[j].weight;
	  
	  double temp_err=temp_weight*sqrt(temp_err1*temp_err1+
					   temp_err2*temp_err2);

	  double reweight =
	    MC_reweight(MC_jet1[j], MC_jet2[j], correlation, nGenerate);
	  
	  double temp_weight_reweight=
	    reweight*rescale_reweight;

	  double temp_err_reweight=temp_weight_reweight*
	    sqrt(temp_err1*temp_err1+
		 temp_err2*temp_err2);

	  DataSet_tau21_2jet_KDE.Fill(make_list<double>
				      (MC_jet1[j].value[0], MC_jet2[j].value[0]),
				      w*temp_weight);

	  DataSet_tau21_2jet_KDE_corr.Fill(make_list<double>
					   (MC_jet1[j].value[0], MC_jet2[j].value[0]),
					   w*temp_weight_reweight);


	  hist["T21_MC"]->Fill(temp_T21, w*temp_weight);
	  hist["T21_MCerr"]->Fill
	    (temp_T21, w*temp_err);

	  hist["T21_MC_corr"]->Fill(temp_T21, w*temp_weight_reweight);
	  hist["T21_MCerr_corr"]->Fill(temp_T21, w*temp_err_reweight);
	}


      //first get the correlation for mpt
      MC_corr=
	MC_mpt_2jet_dR.GenerateFull(make_list<double>(dR));
      
      correlation=MC_correlation(MC_corr);

      //different correlation from 2jet
      MC_corr=
	MC_mpt_2jet_pt1_dR.GenerateFull(make_list<double>(log(pt[0]/50.), dR));

      double correlation_pt1=MC_correlation(MC_corr);
      
      
      //now generate events for mpt
      vector<PDF::MC_Event> MC_mpt1=
	MC_mpt_pt.GenerateGrad(make_list<double>(log(pt[0]/50.)), nGenerate);

      vector<PDF::MC_Event> MC_mpt2=
	MC_mpt_pt.GenerateGrad(make_list<double>(log(pt[1]/50.)), nGenerate);


      //compute reweights
      rescale=0;
      rescale_reweight=0;
      double rescale_reweight_pt1=0;

      rescale=MC_rescale(MC_mpt1, MC_mpt2, correlation, rescale_reweight);
      MC_rescale(MC_mpt1, MC_mpt2, correlation_pt1, rescale_reweight_pt1);
      
      for(unsigned int j=0; j<MC_mpt1.size(); ++j)
	{
	  double mass1=exp10(-MC_mpt1[j].value[0])*pt[0];
	  double mass2=exp10(-MC_mpt2[j].value[0])*pt[1];

	  hist["m1_MC"]->Fill(mass1, w*MC_mpt1[j].weight);
	  hist["m1_MCerr"]->Fill(mass1, w*MC_mpt1[j].err);
	  hist["m2_MC"]->Fill(mass2, w*MC_mpt2[j].weight);	
	  hist["m2_MCerr"]->Fill(mass2, w*MC_mpt2[j].err);
	}

      //fill MJ now
      for(unsigned int j=0; j<MC_mpt1.size(); ++j)
	{
	  double temp_weight=MC_mpt1[j].weight* MC_mpt2[j].weight*rescale;
	  double temp_err1=0;
	  double temp_err2=0;

	  double reweight =
	    MC_reweight(MC_mpt1[j], MC_mpt2[j], correlation, nGenerate);

	  double reweight_pt1 =
	    MC_reweight(MC_mpt1[j], MC_mpt2[j], correlation_pt1, nGenerate);

	  double temp_weight_reweight=
	    reweight*rescale_reweight;

	  double temp_weight_reweight_pt1=
	    reweight_pt1*rescale_reweight_pt1;

	  if(temp_weight>0)
	    {
	      temp_err1=MC_mpt1[j].err/MC_mpt1[j].weight;
	      temp_err2=MC_mpt2[j].err/MC_mpt2[j].weight;
	    }

	  double temp_err=temp_weight*sqrt(temp_err1*temp_err1+
					   temp_err2*temp_err2);

	  double temp_err_reweight=temp_weight_reweight*
	    sqrt(temp_err1*temp_err1+
		 temp_err2*temp_err2);

	  double temp_err_reweight_pt1=temp_weight_reweight_pt1*
	    sqrt(temp_err1*temp_err1+
		 temp_err2*temp_err2);
	  
	  double mass1=exp10(-MC_mpt1[j].value[0])*pt[0];
	  double mass2=exp10(-MC_mpt2[j].value[0])*pt[1];

	  DataSet_mpt_2jet_KDE.Fill(make_list<double>
				    (MC_mpt1[j].value[0], MC_mpt2[j].value[0]),
				    w*temp_weight);

	  DataSet_mpt_2jet_KDE_corr.Fill(make_list<double>
					 (MC_mpt1[j].value[0], MC_mpt2[j].value[0]),
					 w*temp_weight_reweight);

	  total_weight+=w;
	  
	  hist["MJ_MC"]->Fill(mass1+mass2, w*temp_weight);
	  hist["MJ_MCerr"]->Fill(mass1+mass2, w*temp_err);

	  hist["MJ_MC_corr"]->Fill(mass1+mass2, w*temp_weight_reweight);
	  hist["MJ_MCerr_corr"]->Fill(mass1+mass2, w*temp_err_reweight);

	  hist["MJ_2jet_MC_corr"]->Fill(mass1+mass2, w*temp_weight_reweight_pt1);
	  hist["MJ_2jet_MCerr_corr"]->Fill(mass1+mass2, w*temp_err_reweight_pt1);

	}


      //now generate toy MC from dijet distributions
      vector<PDF::MC_Event> MC_mpt12=
	MC_mpt12_pt12.GenerateFast
	(make_list<double>(log(pt[0]/50.), log(pt[1]/50.)), nGenerate);

      
      //loop over events
      for(unsigned int j=0; j<MC_mpt12.size(); ++j)
	{
	  double mass1=exp10(-MC_mpt12[j].value[0])*pt[0];
	  double mass2=exp10(-MC_mpt12[j].value[1])*pt[1];

	  hist["MJ_2jet_MC"]->Fill(mass1+mass2, w*MC_mpt12[j].weight);
	  hist["MJ_2jet_MCerr"]->Fill(mass1+mass2, w*MC_mpt12[j].err);
	}
      
      
      //generate 4jet jet mass
      nGenerate=500;

      MC_jet1=
	MC_mpt_pt.GenerateFast(make_list<double>(log(pt[0]/50.)), nGenerate);

      MC_jet2=
	MC_mpt_pt.GenerateFast(make_list<double>(log(pt[1]/50.)), nGenerate);
      
      vector<PDF::MC_Event> MC_jet3=
	MC_mpt_pt.GenerateFast(make_list<double>(log(pt[2]/50.)), nGenerate);
      
      vector<PDF::MC_Event> MC_jet4=
	MC_mpt_pt.GenerateFast(make_list<double>(log(pt[3]/50.)), nGenerate);


      //get the rescale factor
      double sum_weight=0;
      for(unsigned int j=0; j<MC_jet1.size(); ++j)
	{
	  sum_weight+= MC_jet1[j].weight*MC_jet2[j].weight*
	    MC_jet3[j].weight*MC_jet4[j].weight;
	}
      
      assert(sum_weight > 0);
      sum_weight=1.0/sum_weight;
      
      for(unsigned int j=0; j<MC_jet1.size(); ++j)
	{
	  double temp_mass=pt[0]*exp10(-MC_jet1[j].value[0]) +
	    pt[1]*exp10(-MC_jet2[j].value[0]) +
	    pt[2]*exp10(-MC_jet3[j].value[0]) +
	    pt[3]*exp10(-MC_jet4[j].value[0]);
	  
	  double temp_weight= MC_jet1[j].weight*MC_jet2[j].weight*
	    MC_jet3[j].weight*MC_jet4[j].weight*sum_weight;

	  if(temp_weight==0) continue;

	  double temp_err1=MC_jet1[j].err/MC_jet1[j].weight;
	  double temp_err2=MC_jet2[j].err/MC_jet2[j].weight;
	  double temp_err3=MC_jet3[j].err/MC_jet3[j].weight;
	  double temp_err4=MC_jet4[j].err/MC_jet4[j].weight;
	  
	  double temp_err=temp_weight*sqrt(temp_err1*temp_err1+
					   temp_err2*temp_err2+
					   temp_err3*temp_err3+
					   temp_err4*temp_err4);


	  hist["MJ_4jet_MC"]->Fill(temp_mass, w*temp_weight);
	  hist["MJ_4jet_MCerr"]->Fill(temp_mass, w*temp_err);
	  
	}
    }
  //end event loop

  cout<<endl;


  //generate MC for correlation studies
  //first grab some more histograms
  hist["t21_corr"]=new TProfile("t21_corr", "#tau_{21} Correlation", 100, 1.0, 2*PI);

  

  for(float i=1.2; i<2*PI; i+= ((2*PI-1.2)/100.))
    {
      vector<PDF::MC_Event> MC_t21_evt=
	MC_tau21_2jet_dR.GenerateFull(make_list<double>(i));

      double correlation_norm=
	MC_correlation_norm(MC_t21_evt);
      
      hist["t21_corr"]->Fill
	(i, correlation_norm );
    }

  total_weight=1.0/total_weight;

  //2D Analysis

  //mpt 2D plot
  //now get the mpt1 vs mpt2 PDF plot
  //get the actual PDF from kernel smoothing
  TProfile2D actual_mpt12("actual_mpt12_profile", "Actual m/p_{T} from Sherpa",
			  40, 0, 2.0, 40, 0, 2.0);

  TProfile2D actual_mpt12_err("actual_mpt12_profile_err",
			      "Actual m/p_{T} from Sherpa (error)",
			      40, 0, 2.0, 40, 0, 2.0);

  TProfile2D KDE_mpt12("KDE_mpt12_profile", "m/p_{T} from KDE",
		       40, 0, 2.0, 40, 0, 2.0);

  TProfile2D KDE_mpt12_diff("KDE_mpt12_profile_diff", "m/p_{T} from KDE",
			    40, 0, 2.0, 40, 0, 2.0);

  TProfile2D KDE_mpt12_corr("KDE_mpt12_profile_corr", "m/p_{T} from KDE",
			    40, 0, 2.0, 40, 0, 2.0);

  PDF PDF_mpt_2jet= DataSet_mpt_2jet.ComputePDF(true);
  
  PDF::MC_Generator MC_mpt_2jet= PDF_mpt_2jet.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT));

  vector<PDF::MC_Event> MC_full=
    MC_mpt_2jet.GenerateFull(make_list<double>());
  
  
  for(unsigned int i=0; i<MC_full.size(); ++i)
    {
      actual_mpt12.Fill(MC_full[i].value[0], MC_full[i].value[1],
			MC_full[i].weight);

      actual_mpt12_err.Fill(MC_full[i].value[0], MC_full[i].value[1],
			    MC_full[i].err);
    }

  //from KDE
  PDF PDF_mpt_2jet_KDE= DataSet_mpt_2jet_KDE.ComputePDF(true);
  PDF PDF_mpt_2jet_KDE_corr= DataSet_mpt_2jet_KDE_corr.ComputePDF(true);
  
  PDF::MC_Generator MC_mpt_2jet_KDE= PDF_mpt_2jet_KDE.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT));

  PDF::MC_Generator MC_mpt_2jet_KDE_corr= PDF_mpt_2jet_KDE_corr.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT));

  MC_full=
    MC_mpt_2jet_KDE.GenerateFull(make_list<double>());
  
  
  for(unsigned int i=0; i<MC_full.size(); ++i)
    {
      KDE_mpt12.Fill(MC_full[i].value[0], MC_full[i].value[1],
		     MC_full[i].weight);

      KDE_mpt12_diff.Fill(MC_full[i].value[0], MC_full[i].value[1],
			  MC_full[i].weight);
    }

  MC_full=
    MC_mpt_2jet_KDE_corr.GenerateFull(make_list<double>());
  
  
  for(unsigned int i=0; i<MC_full.size(); ++i)
    {
      KDE_mpt12_corr.Fill(MC_full[i].value[0], MC_full[i].value[1],
			  MC_full[i].weight);
    }


  //now get the TGraph
  KDE_mpt12.Scale(100*100);
  KDE_mpt12_diff.Scale(100*100);
  KDE_mpt12_corr.Scale(100*100);
  actual_mpt12.Scale(100*100);
  actual_mpt12_err.Scale(100*100);

  KDE_mpt12_diff.Add(&actual_mpt12, -1);

  KDE_mpt12.Write();
  KDE_mpt12_diff.Write();
  KDE_mpt12_corr.Write();
  actual_mpt12.Write();
  actual_mpt12_err.Write();

  TGraph2D KDE_mpt12_graph(&KDE_mpt12);
  TGraph2D KDE_mpt12_graph_corr(&KDE_mpt12_corr);
  TGraph2D actual_mpt12_graph(&actual_mpt12);
  
  KDE_mpt12_graph.SetName("KDE_mpt12");
  KDE_mpt12_graph_corr.SetName("KDE_mpt12_corr");
  actual_mpt12_graph.SetName("actual_mpt12");

  
  KDE_mpt12_graph.Write();
  KDE_mpt12_graph_corr.Write();
  actual_mpt12_graph.Write();  


  //tau21 2D plot
  //now get the mpt1 vs mpt2 PDF plot
  //get the actual PDF from kernel smoothing
  PDF PDF_tau21_2jet= DataSet_tau21_2jet.ComputePDF(true);

  TProfile2D actual_tau21("actual_tau21_profile", "Actual #tau_{21} from Sherpa",
			  40, 0, 1.0, 40, 0, 1.0);

  TProfile2D actual_tau21_err("actual_tau21_profile_err", 
			      "Actual #tau_{21} from Sherpa (error)",
			      40, 0, 1.0, 40, 0, 1.0);

  TProfile2D KDE_tau21("KDE_tau21_profile", "#tau_{21} from KDE",
		       40, 0, 1.0, 40, 0, 1.0);

  TProfile2D KDE_tau21_diff("KDE_tau21_profile_diff", 
			    "#tau_{21} from KDE (difference)",
			    40, 0, 1.0, 40, 0, 1.0);

  TProfile2D KDE_tau21_corr("KDE_tau21_profile_corr", "#tau_{21} from KDE",
			    40, 0, 1.0, 40, 0, 1.0);

  
  PDF::MC_Generator MC_tau21_2jet= PDF_tau21_2jet.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT));

  MC_full=
    MC_tau21_2jet.GenerateFull(make_list<double>());
  
  
  for(unsigned int i=0; i<MC_full.size(); ++i)
    {
      actual_tau21.Fill(MC_full[i].value[0], MC_full[i].value[1],
			MC_full[i].weight);

      actual_tau21_err.Fill(MC_full[i].value[0], MC_full[i].value[1],
			    MC_full[i].err);
    }

  //from KDE
  PDF PDF_tau21_2jet_KDE= DataSet_tau21_2jet_KDE.ComputePDF(true);
  PDF PDF_tau21_2jet_KDE_corr= DataSet_tau21_2jet_KDE_corr.ComputePDF(true);
  
  PDF::MC_Generator MC_tau21_2jet_KDE= PDF_tau21_2jet_KDE.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT));

  PDF::MC_Generator MC_tau21_2jet_KDE_corr= PDF_tau21_2jet_KDE_corr.Generator
    (make_list<PDF::MC_Flag>
     (PDF::OUTPUT, PDF::OUTPUT));

  MC_full=
    MC_tau21_2jet_KDE.GenerateFull(make_list<double>());

  for(unsigned int i=0; i<MC_full.size(); ++i)
    {
      KDE_tau21.Fill(MC_full[i].value[0], MC_full[i].value[1],
		     MC_full[i].weight);

      KDE_tau21_diff.Fill(MC_full[i].value[0], MC_full[i].value[1],
			  MC_full[i].weight);
    }

  MC_full=
    MC_tau21_2jet_KDE_corr.GenerateFull(make_list<double>());

  for(unsigned int i=0; i<MC_full.size(); ++i)
    {
      KDE_tau21_corr.Fill(MC_full[i].value[0], MC_full[i].value[1],
			  MC_full[i].weight);
    }


  //now get the TGraph
  KDE_tau21.Scale(100*100);
  KDE_tau21_diff.Scale(100*100);
  KDE_tau21_corr.Scale(100*100);
  actual_tau21.Scale(100*100);
  actual_tau21_err.Scale(100*100);

  KDE_tau21.Write();
  KDE_tau21_diff.Add(&actual_tau21,-1);
  KDE_tau21_diff.Write();
  KDE_tau21_corr.Write();
  actual_tau21.Write();
  actual_tau21_err.Write();

  TGraph2D KDE_tau21_graph(&KDE_tau21);
  TGraph2D KDE_tau21_graph_corr(&KDE_tau21_corr);
  TGraph2D actual_tau21_graph(&actual_tau21);
  

  KDE_tau21_graph.SetName("KDE_tau21");
  KDE_tau21_graph_corr.SetName("KDE_tau21_corr");
  actual_tau21_graph.SetName("actual_tau21");

  
  KDE_tau21_graph.Write();
  KDE_tau21_graph_corr.Write();
  actual_tau21_graph.Write();  


  //cout<<"weights thrown away: "<<low_weight/total_weight*100<<" %"<<endl;
  
  //write all histograms
  for(it=hist.begin(); it!=hist.end(); ++it)
    {
      //see if the name has MC but no MCerr
      if(it->first.find("MC")!=string::npos &&
	 it->first.find("MCerr")==string::npos)
	{
	  //try to see if the corresponding MCerr histogram exists
	  string str=it->first;
	  str.replace(str.find("MC"), 2, "MCerr");
	  TH1* ch=hist[str];
	  if(ch!=NULL)
	    {
	      //found a corresponding histogram
	      for(int i=1; i<=it->second->GetNbinsX(); ++i)
		{
		  it->second->SetBinError
		    (i, ch->GetBinContent(i));
		}

	    }
	}
      
      it->second->Write();
    }

  outfile.Close();
  cout<<"Program ends."<<endl;
  return 1;
}
  
